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The power of matrix product states to describe infinite-size translational-invariant critical spin 
chains is investigated. At criticality, the accuracy with which they describe ground state properties 
of a system is limited by the size x of the matrices that form the approximation. This limitation is 
quantified in terms of the scaling of the half-chain entanglement entropy. In the case of the quantum 
Ising model, we find ~ | log x with high precision. This result can be understood as the emergence 
of an effective finite correlation length ruling all the scaling properties in the system. We produce 
six extra pieces of evidence for this finite-x scaling, namely, the scaling of the correlation length, 
the scaling of magnetization, the shift of the critical point, the scaling of the entanglement entropy 
for a finite block of spins, the existence of scaling functions and the agreement with analogous 
classical results. All our computations are consistent with a scaling relation of the form ~ X") 
with K = 2 for the Ising model. In the case of the Heisenberg model, we find similar results with 
the value k ~ 1.37. We also show how finite-x scaling allows to extract critical exponents. These 
results are obtained using the infinite time evolved block decimation algorithm which works in the 
thermodynamical limit and are verified to agree with density matrix renormalization group results 
and their classical analogue obtained with the corner transfer matrix renormalization group. 



I. INTRODUCTION 

The exact solution of the dynamics of quantum phys- 
ical systems is often too hard or impossible to compute. 
It is then necessary to resort to approximation schemes 
and numerical simulations, as in the case of QCD, the 
theory of strong interactions, to gain some insight into 
the physics of the theory under study. These numerical 
simulations are implemented using some clever algorithm 
that exploits the understanding of the quantum interac- 
tions at work. It may then be difficult to separate what 
is the absolute limitation inherent to the nature of the 
approximation from what is an artifact of the specific 
algorithm employed. 

We can elaborate further this idea in the case of one- 
dimensional translational invariant spin chains. There, 
the recent algorithms introduced by Vidal based on the 
explicit use of Schmidt decompositions 1] have been 
shown to deliver identical results to the very successful 
Density Matrix Renormahzation Group (DMRG)[1,|1,|1|- 
Actually, these two apparently wide-apart algorithms 
agree because they come down to represent the coeffi- 
cients of a quantum state as a product of matrices, that 
is a Matrix Product State (MPS) [1, i, 0] 

I*) = ^ ti[Ai{si)...AN{sN)]\si...SN), (1) 

S1...SJV 

where Si labels a basis for the local degree of freedom 
('spin') of particle i, and where the ylj(si)'s are matrices 
of some fixed finite size, Xj a-nd N is the number of sites 
in the chain which will be taken to be infinite [33] . Under 
the assumption that the above mentioned algorithms do 
find a faithful description of the sought state, consistent 
with the MPS structure, we can forget about their de- 



tails and describe their results as a consequence of the 
properties of MPS states. 

In this paper, we shall investigate what is the limita- 
tion attached to the use of the MPS approximation for 
infinite one-dimensional translational invariant quantum 
systems. It is important to note that we address infinite 
systems in order to avoid the presence of any finite size 
effect. Consequently, any departure of MPS results from 
the exact ones is expected to be due to the very nature 
of the MPS representation and must necessarily relate to 
the finite matrix-size x that can be handled in practice. 

Let us now present a brief summary of our main re- 
sult. We need first to recall the basic construction of 
the Schmidt decomposition for any state in a bipartite 
Hilbert space H = Ha ^Hb, 

nun{dim'HA,'HB) 

l*)= E A„|^f)|0, (2) 

where EJA.!' = 1, = (v'fl^^) = S^p. The 

amount of entanglement (quantum correlations) between 
A and B can be quantified in terms of the von Neumann 
entropy of part A (or B): 

^(PA) = -E^'logA2 =5(pb) . (3) 

This entanglement entropy in an infinite chain is known 
to obey scaling properties. At a critical point [34], the 
entanglement of a block of size L with the rest of the 
chain scales as Q 

SiL)^^logL, (4) 

where c is the central charge associated with the univer- 
sality class of the quantum phase transition. In particu- 
lar, we can take party A to be the left half of the chain 
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with L = N/2 sites and party B to be the right half with 
the remaining sites. It is clear that the entanglement of 
half of the chain with the other half will diverge as N 
goes to infinity. More precisely, if we consider a system 
with open boundary conditions, the following diverging 
behavior is expected 

^(infinite half — chain) ^^^5° ^ log — . (5) 

Asymptotically, for very long chains, the half chain en- 
tropy is only half of the entropy of a block with the 
same size. This can be understood by noticing that the 
block has two boundaries available to establish correla- 
tions with the rest of the chain, whereas a half chain 
only has one. We may now wonder how much of this 
infinite amount of entanglement is captured by the MPS 
approximation. For a system in an MPS with matrices 
of size X: S{pa) is trivially bounded by logx- It is thus 
obvious that an MPS with matrices of finite size cannot 
describe exactly the behavior of an infinite system at the 
critical point but we may try to find the exact amount of 
entanglement which is captured. 

We have found that the quantitative entanglement sup- 
port of MPS at criticality obeys the following scaling law 
for the quantum Ising model 

•^x^^logX (6) 

with a remarkably high precision. 

This effective saturation of the entropy can be under- 
stood in an elegant way as the emergence of a finite corre- 
lation length a fact that was first analyzed in Ref. @ 
in the context of DMRG calculations for gapless systems. 
To complete the connection we use the known result 
that, near criticality, entanglement entropy is expected 
to be saturated by S* ~ |log^. Typical values of the 
central charge are c = 1/2 for the Ising model and c = 1 
for the Heisenberg model. 

Thus, our result hints at the Gnite-x scaling relation 

= with K = 2, (7) 

for the quantum Ising model. Moreover, we shall find 
this relation to be fully consistent with many other scal- 
ing properties in the system. In some sense we may argue 
that the finite matrix-size x inherent to the MPS approx- 
imation works as a probe of the universality class of the 
quantum phase transition which is investigated, a fact 
which is analogous to the well-known finite-size scaling 
for finite systems [TT| . 

Our results will be mainly numerical obtained with a 
specific technique. The best MPS approximation to a 
given state can be obtained using different algorithms, 
DMRG being the most popular choice. Nevertheless, the 
recentl y pr oposed infinite time-evolving block decimation 
iTEBD [Tz, [ij] turns out to be particularly suited to ad- 
dress infinite quantum systems. This algorithm exploits 



translational invariance, makes the programming quite 
simple and, for our purposes, runs faster than the com- 
monly used finite size DMRG. Yet, we have verified that 
the results we are presenting here can be obtained using 
DMRG. We are therefore led to believe that our find- 
ings are intrinsic to the MPS representation and are not 
really sensitive on the precise algorithm used to get an 
approximation of the ground state. 

We also have compared our results with the corre- 
sponding classical ones when available. The agreement 
emerging from this comparison is a hint that the scal- 
ing properties we are facing of MPS could be a general 
phenomenon for quantum phase transitions studied with 
tensor network techniques (of which the MPS is just a 
possible choice). 

We would like to stress that our goal is to settle the 
scaling properties inherent to the MPS approximation. 
For that purpose, we do not need to work with MPS's 
with matrices of very large size x we reach the 

scaling region. This region, for the case we study is de- 
fined by 

» « (8) 

where a = 1 is the lattice spacing. Hence, depending on 
the value of n, the scaling region can be attained with 
very modest values of x- 

The paper is organized as follows. In section |TT] we 
discuss the origin of a finite-x scaling relation. Then, 
in section IIIII we collect numerical evidence supporting 
its validity. In section IIVI we show that a similar scal- 
ing relation is expected for the Heisenberg model. Some 
applications of finite-x scaling are briefly discussed in 
section |Vl Namely, we will show how to extract criti- 
cal exponents from finite-x scaling. We summarize our 
results in section IVIl Details regarding the iTEBD algo- 
rithm, its convergence and some improvements we have 
implemented are presented in the Appendices. 



II. FINITE X SCALING 

Phase transitions are usually detected through a lo- 
cal order parameter that discriminates between the two 
phases separated by the critical point. Let us consider a 
concrete example, the infinite quantum Ising model in a 
transverse field [ij] 

■i 

The phase transition of this model is driven by the trans- 
verse magnetic field, A. The x-magnetization plays the 
role of an order parameter and scales as M = (af) ^ 
\\2 _ A*2|i/8 near the critical point A* = 1 [15]. 

We expect that, at criticality, a description of the 
ground state of H in terms of a finite x MPS blurs a phase 
transition smooth. For instance a diverging correlation 
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length at A* = 1 is replaced by a peak for the value of 
at some value A* of the transverse field (A* = A*^^^). 
Indeed the correlation length of an MPS is usually finite 

[i,[ii[33. 

The value of the peak, and its position A* should 
be dictated by a scaling relation of the following type 



(10) 



Let us briefly argue why this should be the case, by show- 
ing how the arguments in Ref . [ll| . formulated for finite 
size scaling, can be adapted to the case of finite x scal- 
ing. If Eq. (jlOp holds, in analogy with what is observed 
in finite systems, the MPS finite x smooths all the diver- 
gences that we would observe in infinite systems at the 
phase transition. They should be transformed to some fi- 
nite anomaly at a % dependent pseudo-critical point A* . 
To see this, we start by noticing that, asymptotically, the 
correlation length depends only on the distance from the 
transition through the universal critical exponent v. 



(11) 



where t = |A — A*|/A*. By reading this relation in the 
opposite direction we gain some further understanding 



(12) 



Given that x cannot be taken to infinity, we are keep- 
ing the system away from criticality. The transition is 
actually shifted to a pseudo phase transition located at 
a different value of the magnetic field A*. There, the 
correlation length does not diverge. By substituting Eq. 

into Eq. we obtain a prediction on how the 

pseudo-critical point should approach the true critical 
point when varying x- 



|A* - A* 
A^ 



X 



(13) 



For a given x, we obtain the effective distance from 
criticality when the system is at its critical point. We can 
hence stick there , at A*, and fix our attention on how 
universal quantities should vary as we change x- We may 
now envisage three different scenarios. When a universal 
quantity diverges approaching the critical point with 
an exponent lo this translates to a divergence at A* in 
term of x ^"^'^ 



Fu{\*) - X 



(14) 



In the case where the universal quantity vanishes when 
approaching the critical point with a given exponent u, 
as is the case for the order parameter, then we should 
have 



F„(A* 



X 



(15) 



As a last case, we consider the possibility of a logarithmic 
divergence, as is the case for the half chain entropy. Then, 



F„(A*) ^Klog(x). 



(16) 



Now we can look for deviations from the critical point. 
Once we have isolated the anomalous contributions to the 
universal quantities we are left with a regular part that, 
if correctly interpreted, does not depend on the size of 
the matrices. In analogy to the finite size case, we call 
this contribution the scaling function for that particular 
universal quantity. An intuitive picture of its origin can 
be obtained by considering again Eq. (fTTj). We consider 
the variable 



(17) 



that, in an infinite system, stays of order one in all the 
critical region, including the phase transition point as 
guaranteed by Eq. (jlip . Away from the critical region, 
where the correlation length attains a finite value, it in- 
creases monotonically with t. When passing to finite x 
systems we break the relation Eq. pT|) . Expressing the 
correlation function in term of x by using Eq. (jlO[) we 
get 



tx 



k/i/ 



(18) 



Values for this variable close to zero, are due to finite 
X effects and can easily be obtained by getting closer and 
closer to the critical point at fixed x- This is the vari- 
able that really quantifies the distance from an infinite 
system. Systems with different x ^-t different t but with 
the same x are indeed at the same distance from the cor- 
responding infinite system. The variable x can thus be 
used to unmask x independent effects induced by forcing 
the system away from its critical behavior. In order to 
do this, however, one should keep in mind that systems 
with different x have also different anomaly strengths as 
described by equations (fT4|) . (fTS]) and (fT6|) . In order to 
unmask x independent effects we should therefore nor- 
malize the results obtained with system with different 
X with their anomalous contributions at the transition. 
For the cases considered in Eq. p^ . and (|16p the 
scaling functions are extracted respectively as 



fu{x) ~ X Fu{x), 



Iu{x) ~ x^" Fu{.x), 



fu{x) 



Fujx) 
Klog(x)' 



(19) 



(20) 



(21) 



We now provide numerical support to the finite-x scal- 
ing. 



III. EVIDENCE FOR FINITE-x SCALING FOR 
THE QUANTUM ISING CHAIN 

The general discussion on finite-x scaling should be 
verified on concrete examples. We present in this section 



the results for the quantum Ising chain in a transverse 
magnetic field in Eq. All our results are obtained 

using the iTEBD algorithm. Some aspects of this tech- 
nique are discussed in the Appendix. 



A. Half-chain entropy 

We first compute the von Neumann entropy for half 
the infinite chain. As mentioned previously, this mea- 
sure of entanglement should diverge with the size of the 
system. Such a divergence cannot be accommodated by 
a finite-x MPS ansatz. Entanglement must circulate via 
the ancillary indices of the matrices that build the ap- 
proximation. For matrices of size Xi entanglement is 
bounded to only span a space of dimension 2^ as ex- 
plained in Ref. ^17], rather than the actual diverging 2^ 
dimensions. Moreover, the eigenvalues in the Schmidt 
decompositions obey some decay law (an exponential de- 
cay, up to degeneracies, is expected from confornial field 
theory), that further decreases the amount of entangle- 
ment that the approximation should support. 

Numerical results for the entanglement entropy for the 
half-chain at A = 1 is shown in Fig. [1] where we have 
plotted as a function of % and found an accurate fit 
to the scaling law 

^x^^logX. (22) 

The remarkable precision of the fit should emerge from 
the absence of constant and y^— corrections. This effect 
was observed in the context of block entropies in Ref. |18j 
and shown absent in other measures of quantum correla- 
tions like the single copy entanglement which is based on 
the largest eigenvalue of the reduced density matrix of a 
subsystem. Conformal symmetry orchestrates a cancella- 
tion of subleading terms coming from all the eigenvalues 
of that reduced density matrix. In the present case, it is 
unclear why corrections are absent in the computation of 
the half-chain entropy at the point A = 1. 



We can now match this scaling to the result of Ref. fld\ 
that states that, away from criticality, log^. Then, 

the hypothesis of finite-x scaling suggests the half-chain 
entropy to behave as 

5x~^logx^ (23) 

In the case of the Ising model, we find 

^^X", with K ~ 2.011(2), (24) 

where we have used that the central charge c is equal 
to 1/2 for the Ising model. The error in our result only 
reflects the quality of the fit. This depends on the use of 
small values of x, where scaling may not be present, and 
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FIG. 1: Entropy as a function of logx at A = 1. 

on possible violations of that scaling. The uncertainty is 
then not representing a faithful systematic error but just 
the order of magnitude of the freedom in the fit. Our goal 
in this paper remains to collect a first consistent estimate 
for what is the actual value of k. 

In practical terms, this result shows that numerical 
exploration of the critical properties should be well de- 
scribed using relatively small MPS. A value of x 20 
describes faithfully correlations up to 400 sites. 

We now consolidate this result by checking its consis- 
tency with the computation of other observables. 



B. Shift of the critical point 

In the vicinity of the critical value A* = 1, the entan- 
glement entropy of half of the Ising chain diverges and the 
magnetization abruptly drops to zero. The best MPS ap- 
proximation to this scenario manages to produce a peak 
in the entropy and sudden drop of the magnetization for 
values of A which are shifted from the infinite chain crit- 
ical value. We label A* g the coupling where the entropy 
presents a peak and the A* J^J the one where the mag- 
netization vanishes abruptly. As in finite size simulation 
schemes Tl] , we have found that both A* 5 7^ A* and 
A* 7^ A*. But we have found that, within the accu- 
racy of our simulation. A* = A* 5 = A* . This can be 
understood as a check of the consistent representation of 
criticality that MPS develop. 

Our results are shown in Fig. [2] for the entropy and 
the magnetization respectively. We can see that (i) the 
amplitude of the shift A* — A* reduces when we increase 
X, (ii) the peak of the entropy rises with increasing Xj 
and that (iii) far from the critical point, modest values 
of X are sufficient to get faithful approximation of the 
ground state (in the sense that the curves obtained for 
different values of x tend to collapse). 

We have checked that the shift of the critical point 
obeys the law (fT5|) . The results are plotted in Fig. [31 As 



expected, the way A* approaches A* is correctly described 



by a power law. Using — 1 we extract: 



2.1(1) 



(25) 



where, again, the error is only reflecting the precision of 
the fit. 

This value is compatible with the value extracted us- 
ing the entanglement entropy. We see, however, that 
this estimation is less precise. This fact is related to the 
difSculties encountered in the determination of A*. In 
principle the finer the scan, the more precise the value of 
A* . However the sharpness of the scan is limited by the 
numerical precision with which we obtain the entropy. 
At some point, entropies of chains with close but differ- 
ent values of A are compatible within their error bars. 
Then, we cannot further refine our scan and should ac- 
cept the obtained precision as the best we can achieve for 
the location of the transition. 
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FIG. 3: Effective critical point as a function of x- As 
discussed in the appendices, the values for X = 2 and 4 were 



obtained with e 



10" 



while for x 



16 and 32 with 
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The errors bar are due to the finite value of s. 



C. Magnetization 
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The drop of the magnetization near the critical point 
obeys scaling laws as discussed previously. We actually 
expect the magnetization at finite x to behave as (A = 

^*) ^ with the Ising critical exponents /3 = 1/8 and 

— 1. We may now take our numerical results and fit k 
in this expected scaling law. In Fig. IH we have plotted 
A/^ as a function of x for the Ising chain at A = 1. By 
fitting our numerical results with a function of the form 



(see Fig. SJ, we obtain: 
K = 2.03(2). 



(26) 



This value of k is in agreement with our two previous 
determinations . 




FIG. 2: Entropy and magnetization around the theoretical 
critical point obtained for x = 2, 4, 8 and 16 using, as de- 
scribed in the appendices, e — 10~^ and after convergence of 
the eighth decimal. The error bars due to the finite value of 
e are smaller than the points size. 
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FIG. 4: Magnetization as a function of x at A = 1. 



D. Block entropy E. Correlation length 



A new consistency check consists in considering the 
entropy of the reduced density operator of a block of L 
contiguous particles. For a critical systems, this entropy 
scales with L as Sl — ^ logi 33- We have observed that, 
for a fixed value of x, this entropy saturates at a distance 
L cf. x'^- We can make a very qualitative assumption on 
the fact that the length at which the entropy saturate 
is of the order of the correlation length and in this way 
use this value as determination of the correlation length. 
It is likely that this qualitative assumption can be made 
rigorous in a renormalization group framework by intro- 
ducing a new scaling field x'^^'^ ■ However this analysis 
is out of the scope of this paper. By using the relation 
(fTO|) with this estimation of the correlation length we can 
have an idea of the magnitude of k. 

In order to compute the entropy of a block of L spins, 
we have used the ideas contained in the work by Ver- 
straete et al. The basic idea is to reconstruct the 

effective new matrix MPS upon successive RG coarse- 
graining transformations. Our results are displayed on 
Table HI where we can see that Sl saturates for L ~ . 
So that we get a further confirmation that, for the Ising 
model, 



K ~ 2.0(1), 



(27) 



in agreement with the previous estimations, though less 
accurate. We observe that for sufficiently large L, Sl 
is approximately equal to two times 5*^ (the half-chain 
von Neumann entropy) calculated at the same A. Let 
us recall that the explanation for this factor 2 is that 
a finite block has two boundaries available to establish 
correlations with the rest of the chain, whereas a half- 
infinite chain only has one. 



All our previous results should be a consequence of the 
emergence of a finite correlation length . This fact was 
first investigated in Ref. 0. We can address this point 
by analyzing the ratio of the two highest eigenvalues of 
the transfer matrix computed from the matrices in 
the MPS. On Fig.®, we have plotted the value of as a 
function of x- To extract the value of the exponent k, we 
have performed a fit to numerical data with a function 
of the type ax'^ with a and k left as free parameters. We 
have found 

K = 2.00(3). (28) 
Again, the consistency of this result with our previous 

dptprminfltinns is ma.m'fpst 
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FIG. 5: Correlation length as a function of the size of the x 
in the case of the Ising model at A = 1. 
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2 


0.2994 


0.4825 


0.5883 


4 


0.3279 


0.5647 


0.6976 


8 


0.3317 


0.6271 


0.7934 


16 


0.3317 


0.6586 


0.8720 


32 


0.3317 


0.6586 


0.9288 


64 


0.3317 


0.6586 


0.9577 


128 


0.3317 


0.6586 


0.9630 


256 


0.3317 


0.6586 


0.9632 


512 


0.3317 


0.6586 


0.9632 


1024 


0.3317 


0.6586 


0.9632 



TABLE I: Entropy of a block of L spins using the ideas con- 
tained in [l3]. We observe that the entropy saturates around 
L ^ X^- Note that the values obtained for the entropy af- 
ter saturation are the double of those obtained for half of the 
chain. This factor of two is due to the fact that here the block 
has two boundaries. 



F. Scaling function for the magnetization 

A further manner to test finite-x scaling is to analyze 
in more detail the magnetization. It follows from the 
scaling analysis in Sect. [n]that depends on x only 
through the product x — x'^^'^t- Therefore, we can plot 
the rescaled magnetization M^{x'^^'^t)x'^^^'^ as a function 
of x'^^'^t for different values of x, assuming the known 
values oi u = 1 and (3 = 1/8 of the Ising universality 
class. In case finite-x scaling is verified, all points should 
lie on the same curve. The quality of this collapse is, 
hence, a function of the correct value of k alone. 

We have scanned k for a broad range of values and 
selected the ones that qualitatively produced a collapse 
of the numerical points onto a single curve. Remarkably, 
we have verified that only for a relatively small interval 
of K values, all the points obtained with this procedure 
lie on the same curve. Whatever small variation outside 
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this interval of k reflects on a sensible spread of the point 
outside the curve. 

Our results are displayed on Fig[6] Again, we find a 
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FIG. 6: Collapse of the rescaled magnetization obtained with 
different MPS using distinct values of k on the scaling func- 
tion: K = 1.9,2.0,2.1 for upper, middle and lower graphs, 
respectively. 



G. Scaling function for the the energy difference 
and comparison with classical results 

In this section, we clarify why studying the deviation 
from the critical point with finite MPS we face a two scale 
problem. This problem is the quantum version of the al- 
ready studied two scale finite-size scaling ansatz in the 
context of classical systems [l^. In our case, the anal- 
ysis is performed considering infinite size systems. The 
first scale is given by the MPS dimension x and the sec- 
ond scale is given by t. We saw in the previous section 
that by correctly treating the two scales, one is able to 
extract universal scaling functions. The universality, im- 
plies that scaling functions however, should not depend 
on which system, among those in the same universal- 
ity class, one decide to consider. This statemet can be 
checked by comparing our results with the one contained 
in Ref. [l9| about the classical Ising model in two di- 
mensions at critical temperature. This system is indeed 
in the same universality class we are considering: the 
two dimensional Ising universality class. The authors 
apply the ideas of the Corner Transfer Matrix (CTM) 
renormalization group [20| to it. This is a technique that 
generalizes the DMRG renormalization ideas and its re- 
lated variational techniques over MPS to a real space 
renormalization algorithm for classical systems. Once the 
CTMRG is applied at critical temperature as in ref. [l^ 
a new scale emerges. This is the inverse of a correlation 
length depending on the dimension of the renormalized 
CTM m. The authors label this scale ^{m). This scale 



exactly corresponds to the scale we are calling here . 
In this way, the authors of ref. , by treating the finite 
size classical system of dimension N studied with a finite 
renormalized CTM of dimension to as a two scale prob- 
lem, extract the value of all critical exponents. Here we 
consider the precise map between our results and the one 
contained there. Following the recipes in Ref. [20| we see 
that, (as already implicit in the identification of scales) 
the classical correspondence of x is the size of the renor- 
malized CTM (called to in both references). We can use 
again the finite size scaling ansatz and map the distance 
from the critical point that we call t to the size of the 
classical system considered in ref. [l^ ( called TV In 
this way the scaling variable Xm = of ref. [19'] is 

related to the one we use x as : Xm = x^ . As for the Ising 
model v ^ \, the results in ref. [l^ should exactly corre- 
spond to ours. We check this claim by consider what in 
our language would be the plot in figure 3 of Ref. [l^ ■ 
It represent the energy difference as a function of x with 
respect to the exact result. We can use the standard map- 
ping between the free energy of a classical system and the 
ground state energy of the corresponding quantum sys- 
tem and compare the plot in figure 3 of ref. [I9j . with our 
results for the ground state energy per bond. To do this 
we plot the difference of the ground state energy with a 
given X and the exact result £;(oo) = -1.27323954. The 
scaling function for the energy difference is obtained by 
plotting {E^{l/x) - Ery:_)x"-I" Again we see that by 
using a value of k = 2 the points obtained with different 
X collapse to a single curve. We also see that the curve 
we obtain has the same shape but a different normaliza- 
tion factor with respect to the one obtained in figure 3 of 
Ref. [l^. This is expected since the normalization fac- 
tors are known to be universal but boundary condition 

□ 




FIG. 7: Collapse of the rescaled energy difference obtained 
with different MPS. The values we use are k = 2, E{oo) = 
-1.27323954. 

We can hence confirm that the scaling observed in this 
work in the case of a quantum phase transition is the ana- 



logue to the one observed in a classical phase transition 
in ref. [11. 



IV. 



EVIDENCE OF FINITE-x SCALING FOR 
THE HEISENBERG CHAIN 



An extensive analysis of the emergence of finite-x scal- 
ing in different models is necessary to gain insight in the 
role of the scaling exponent k. Here, we only make a first 
step and explore the Heisenberg spin 1/2 Hamiltonian 



(29) 



-r 



-r 



-1 1 — I — r- 



o Data 

■ ■ ■ Best Fit with Z=[2:44] ; a=0.28(l), b=0.243(3) 
- ■ Best Fit with X=[I6:44] ; a=0.356(9), b=0.226(2) 



S(x) = a + bLog,(x) 



We may conjecture that n should only vary with the uni- 
versality class of the model considered. To assess the 
new value of k, we consider the scaling of the half chain 
entropy since this strategy provided very precise deter- 
mination in the Ising case. 

We then follow the same steps as described for the Ising 
case and we take the central charge to be c = 1. By fitting 
the numerical data with a curve of the type a -I- 6 log x and 
using the actual value of the central charge we obtain, as 
observed on Fig. [51 



K = 1.36(2). 



(30) 



Let us note that the fit now includes a non-zero intersect. 
This was absent in the Ising case. 

This result can be checked for consistency in a similar 
way as the results presented for the Ising model. Here, 
we present as a further piece of evidence for finte-x scal- 
ing the scaling of the correlation length as computed from 
the ratio of the largest eigenvalues of the transfer matrix. 
As shown in Fig. [9l the numerical data are described cor- 
rectly by a law of the type in Eq. (jlOp with an exponent 



AC = 1.38(2). 



(31) 



FIG. 8: Entropy as a function of log x for the Heisenberg 
model. Data have been fitted with a function of the type 
a + 61og(x) with a and h free parameters. The results of 
the factor b for the fit in the x interval from 16 to 44 is 
h = 0.226(2). 
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Data 

Best Fit with x=[2:44]; a=().43(3), K=1.38(2) 
Best Fit with ^=[16:44]; a=0.45(5), K=1.37(3) 
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FIG. 9: Correlation length as a function of x in the case of the 
Heisenberg model. This behavior can be correctly described 
by a relation of the type 1101 with an exponent k = 1.38(4). 
The fit has been performed in the x interval 20 — 44 



Both determinations in Eq. ([30| and (|3ip are compati- 
ble and support the value k, ~ 1.37(2), which depends on 
the universality class of the model under discussion [37|. 



V. APPLICATIONS OF FINITE-x SCALING 

As in the case of finite size scaling, we can use finite x~ 
scaling to extract critical exponents. The ideal strategy 



is the one that does not rely on the knowledge of the 
position of the finite x pseudo critical point. This is so 
because the determination of the pseudo critical point is 
very delicate. Any small error in it propagates to the 
determination of critical exponents as we explicitly saw 
when dealing with the determination of k. 

Keeping this in mind, we can envisage two different 
scenarios: a first simple scenario, as the one of the Ising 
model, when we know a priori the location of the phase 
transition. In this case, in order to extract the critical 
exponents we proceed as follows; i)Extract the value of 
K, by studying the behavior of at the critical point, 
ii) Extract all the ratios a/v where here a represents a 
generic critical exponent by studying universal quantities 
as function of x at the phase transition using the value 
of K, obtained in i) . iii) Extract the value of v (and hence 
a from the ratio obtained in ii) )by studying the deriva- 



tives of the universal quantities with respect to t. The 
second scenario and by far the most unfavorable and fre- 
quent is the one where we do not know the location of the 
critical point. In this case, we need to adapt some strat- 
egy known from finite size analysis to extract its value 
if we want to apply finite x-scaling to the transition. A 
possibility is obtained by considering the techniques of 
22] (more efficient methods can be found i.e. in Ref. 
23l . I2JI). A review of this method for the case of finite 
size scaling is contained in ref. [1^. We adapt it in the 
following way: we iteratively obtain estimates of k and 
the critical point by considering the behavior of the cor- 
relation length as a function of increasingly big x- Once 
these estimates converge to a fixed value, we can use the 
obtained values for k and for the critical point to repeat 
the steps from i) to iii) of case one. In this way we extract 
all the other critical exponents. The main source of error 
in all these determinations is, as in the case of the ffiiite 
size scaling, the existence of scaling violations that we do 
not analyze in this work. However, even without taking 
the scaling violations into account, we think that the ex- 
tracted exponent should be much more accurate than the 
ones obtained with standard techniques. To justify this 
statement, we review what we mean by standard tech- 
niques for extracting critical exponents with an infinite 
MPS by considering again the case of the Ising model. 

We can extract the value of the exponent ly by study- 
ing the behavior of the correlation length at fixed x when 
we approach the phase transition. We expect that far 
enough from the region where finite-x effects appear, a 
modest value of x should provide a faithful description 
of the Ising ground state. The correlation length should 
obey a law of the type (A — A*)""^. Fitting the data with 
this function and leaving A* and v as free parameters we 
obtain an estimate of both A* (the phase transition point) 
and ly. We also expect that due to systematic errors in- 
duced by the fitting procedure (the difficult point is to 
locate the correct window of A values for which we should 
perform the fit) these estimates would have a slight de- 
pendence on X and should converge to the exact A* and 
1/ for X large enough. In Fig. [10] we show the results of 
such study, again for the Ising model. We extract as best 
estimate of ly in the case of x = 16 



ly ~ 1.00(5). 



(32) 



See also Ref. [19|,|26|,|27[ and references therein to see how 
in the case of finite chains described with MPS, finite size 
scaling can be used as an alternative to extract critical 
exponents . 

A similar strategy can be used to extract the (3 critical 
exponent. Again, working slightly away from criticality, 
the scaling of the magnetization is very nicely fitted with 



P ~ .1250(1) 
as shown in Fig. [TlJ 



(33) 



o Data for X=16 
- ■ Best Fit 
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FIG. 10: Fine tuning of the correlation length around the 
critical point for x = 16 and e = 0.1. Note that the points 
are not equally spaced. Inset: Values obtained for f by fitting 
magnetic field window of different sizes (all starting at A = 
0.90 and finishing at the point in the x-axis). We notice a 
good region of stability which can be used to extract our best 
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FIG. 11: Fine tuning of the magnetization around the critical 
point for x = 16 and e = 0.1. Note that the points are 
not equally spaced. Inset: Values obtained for f3 by fitting 
magnetic field window of different sizes all starting at A = 
0.90. We clearly notice a good region of stability which can 
be used to extract our best estimate for the /3 exponent. 



In addition to the exponent /3, one can consider the 
exponent r\ by studying the behavior of the two point 
correlation function of the order parameter a^c- Both ex- 
ponent are related via the hyper scaling relation: /3 — 
{d — 2 + rj)/2, where in this case d = 2 as we are consider- 
ing the universality class of the classical two dimensional 
Ising model. 

This relation implies that \i (3 — 1/8, then rj should 
be 1/4. We checked this for consistency. We plot the 
two point correlation function of the order parameter 
as a function of the distance in Fig. [121 In a log-log 
plot, an algebraic decay such r^^ is seen as a straight 
line. We plot this straight lines together with the cor- 
relations functions obtained for the MPS at the phase 
transition with x = 16,32,64. We appreciate how the 
range for which the correlations reproduce the exact re- 
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suit increases with the matrix dimension. Once the range 
of distances is correctly selected, a fit to a power law in 
the case of correlation function of the x = 64 MPS at 
A = 1 produce the following best estimate for rj 



7] ~ .24800(25). 



(34) 



Again, this result only reflects the quality of the fitting 
strategy. 



X=16 
- ■ X = 32 
■ - X = 64 
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FIG. 12: Study of the order parameter two point correlation 
function at A = 1 for x = 16, 32, 64 and e = 0.01 compared 
with the expected exact behavior r""'^^. We note that the 
range of distances for which there is good agreement between 
the numerical correlation function and the exact result in- 
creases with X as expected. Inset: Results of fits with a power 
law of the type ar~^ for the case of x = 64 in the r windows 
for which the extracted correlation functions agree with the 
analytical results. 



VI. CONCLUSIONS 

The amount of entanglement supported by the MPS 
approximation is limited by the size x of matrices 
that form the ansatz. We have studied numerically this 
issue and found that all observables we have considered 
approach their exact values at criticality obeying scaling 
laws in x- The case of the quantum Ising chain in a trans- 
verse field is consistently described by an effective finite 
correlation length that scales as = x'^j with k ~ 2. 
Most of the results presented here were related to the 
Ising model, but the numerical work we have performed 
shows that our findings are qualitatively valid for other 
models, such as the Heisenberg model where our calcu- 
lations indicate that k ~ 1.36. Interestingly, the value of 
K seems to be model-dependent. 

In the case of the Ising model, it is specially interesting 
to note the accurate fit of the half-chain entropy to 5 ~ 
^ log X at A = 1 with no constant or important subleading 
corrections. This effect is not present for the Heisenberg 
model. 



All our numerical results were found using the iTEBD 
algorithm and checked to agree with standard DMRG 
[a [2^. It would be, in principle, possible to use other 
algorithms as a brute force minimization of energy in the 
space of matrices in the MPS structure. Such an ap- 
proach may fail due to the proliferation of local minima. 
Somehow, DMRG and iTEBD manage to circumnavigate 
local minima an find the absolute minimum within the 
approximation. 

We have also checked that the scaling we encounter 
here coincides with the emergence of a second scale in 
some treatment of classical phase transition as pointed 
out in Ref. [H. 

This correspondence is a hint that this phenomenon is 
quite general and appears whenever one tries to approx- 
imate operators with an infinite rank (such as the CTM 
or the half chain reduced density matrix) with finite rank 
operators Therefore it is likely that scaling is not strictly 
related to the MPS representation of the ground state. 
We are currently investigating this issue by repeating a 
similar study to the one presented here with different 
tensor network representations (29j . 

With the same reasoning, we expect finite-x scaling to 
appear for some generalizations of MPS, such as Tensor 
Product States ^Sfli] also known as Projected Entangled 
Pairs States (sij . It remains an open problem to derive 
the scaling relation analytically for exactly solvable mod- 
els. 
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APPENDIX A: ERROR CONTROL AND 
CONVERGENCE ISSUES WITH THE ITEBD 
ALGORITHM 



In this section, we wish to address the reliability of 
the data output by the iTEBD algorithm. Let us start by 
reminding the main features of this algorithm. A more 
technical presentation can be found in [12.] . 

The iTEBD algorithm aims at finding the ground state 
energy per particle of a Hamiltonian of the form 



H = 



(Al) 



where hi represents a two-spin next-neighbor interaction 
term. This algorithm is based on the following identity. 
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valid for any gapped Hamiltonian: 

1*3) =A/' lim e-^^|vI/o). (A2) 

r — *oo 

That is, a ground state of H can be obtained by evolv- 
ing some initial state in imaginary (Euclidean) time 
whenever H has a gap above the ground state and 
(^ol^'g) 7^ 0. For many Hamiltonians of interest, though, 
Ea. (|A2p cannot be used as such. Rather, one computes 
the following sequence until convergence is attained: 

^-,+1 = £,(e,i/)*,/||£,(e,i/)*,||, (A3) 

where e is some tunable parameter such that £i(e) ~ 
e"*^^ for e small enough. In the iTEBD algorithm, £i(e, H) 
is decomposed into 

£,(e,i/) = ?,?,(£, i?), (A4) 

where the factors appearing in the last expression cor- 
respond each to a different approximation that makes 
numerical computations tractable: 

• i) The first factor ?i(e, H) comes from using a cut 
off Suzuki- Trotter expansion [H in order to approx- 
imate the action of e"*^^ by a product of two-body 
operators. (As a result, the form of 'Ji depends on 
i.) The error introduced by truncating the Suzuki- 
Trotter expansion vanishes when e ^ 0. We call 
this error finite time step error. 

• ii) The second factor CP^ is a projector that ap- 
proximates 3'i(e, i/)^i by an MPS with matrices 
of some prescribed finite size x- This approxima- 
tion is made in order to have an efficient description 
of the state at each step of the sequence (jA3|) . In- 
deed, both storing of ^i+i and the computation of 
the mean value of a local operator now takes a time 
that is polynomial in x UM- This approximation 
boils down to limiting the amount of correlations 
present in the system. We will call truncation error 
the error due to this approximation. 

• iii) The operator e^'^^ is not unitary, and as a result 
'S'i7i{t, H) neither is. This non-unitarity have small 
spurious effects that we can safely neglect [ssj . The 
third operator, Q^, does exactly this job, producing 
what we call an orthonormalization error. 

In order to study the time-step error, we have applied 
the iTEBD algorithm to obtain an MPS approximation of 
the ground state of the quantum Ising chain with matri- 
ces of size X equal to 2. The reason why we have chosen 
to discuss the time-step error with such a small value of 
X is that it is most illustrative. For various values of e 
ranging from 10~^ to 10^^, we have computed the be- 
havior of the ground state energy and the half-chain von 
Neumann entropy. It is natural to test the performance 
of the algorithm looking at the ground state energy since 
it is designed to minimize this quantity. It is less obvious 



A S1-S5 S2-S5 S3-S5 S4-S5 

0.9 9124 921 91 8 

1.0 42495 4203 416 38 

1.071 368647 35489 3501 318 

1.072 69632 6951 689 62 

1.073 69200 6908 684 62 

1.1 58718 5871 581 53 



TABLE II: Convergence of the entropy as a function of e for 
some values of A. The table shows the difference of the entropy 
found using a given e with our best simulations corresponding 
to e = 10^"" ( SI is the value of the half-chain entropy obtained 
for e — 10~^, S2 the values obtained for e = 10~^ and so on). 
All entries in this table should be multiplied by 10~*. 

A AE AS 

0.5 < 0.1 39 

0.7 < 0.1 706 

0.8 < 0.1 2521 

1.0 7 42495 

1.071 41 368647 

1.072 34 69632 

1.073 34 69200 

1.1 29 58718 

1.4 7 14226 

1.5 4 9807 



TABLE III: AE and AS' corresponding to the difference be- 
tween the values obtained for the entropy and energy when 
using £ = 10~^ and e — 10~^ for x = 2. This gives an estima- 
tion of the error due to e when using 10~^ as its value. Note 
that the errors increase around the critical point and that the 
errors in the entropy are much greater than the ones in the 
energy. All entries of this table should be multiplied by 10~*. 



why we also looked at the half-chain entropy. We will 
explain it shortly. 

In principle, the smaller e, the more accurate the de- 
scription of the state. But small values of e also increase 
the number of time steps necessary to guarantee conver- 
gence of the simulation. One way to proceed, in order to 
correctly choose e, is as follows: (i) run a simulation with 
a rather large value of e, ei, and get an estimate of the 
energy and the entropy, (ii) Repeat the simulation with 
a smaller value of e, £2, and compare the resulting en- 
ergy and entropy with those of the simulation at ei . (iii) 
If the results are close enough (according to a predeter- 
mined margin), stop the simulation. Otherwise, repeat 
with smaller values of e until convergence is attained. 

On Table |TT1 we report on the convergence of the von 
Neumann entropy, as a function of e, for various values 
of A, while Table IIIII shows the difference of the values 
for the energy (resp. entropy) for e = 0.1 and e = 10~^. 
We interpret these differences as an estimation of the 
finite time step error at e. (The results for the energy 
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with e = 0.1 and e = 0.01 are already identical up to 
8 decimals. This is why we have not shown them.) We 
observe from these tables that the error on the entropy 
is about ten times larger than that on the energy and 
that both increase around the pscudo critical point A*. 
Simulations with x = 4 and x = 8 show similar results. 
If we now compare the values of the energy and entropy 
yielded by our simulations with the exact values for an 
infinite chain (Table IIVI and Table |V| , we see that the 
errors are larger in the vicinity of the critical point and 
that, as expected, they decrease as we decrease e. 



A 


Exact 


X = 2 


X = 4 


X = 8 


X= 16 


0.5 


1.06354440 


33 


< 0.1 


< 0.1 


< 0.1 


0.6 


1.09223858 


172 


< 0.1 


< 0.1 


< 0.1 


0.7 


1.12682867 


745 


2 


< 0.1 


< 0.1 


0.8 


1.16780951 


2978 


12 


< 0.1 


< 0.1 


0.9 


1.21600091 


12173 


126 


< 0.1 


< 0.1 


1.0 


1.27323954 


69712 


4683 


261 


15 


1.1 


1.34286402 


146576 


1642 


6 


< 0.1 


1.2 


1.41961927 


77696 


416 


< 0.1 


< 0.1 


1.3 


1.50082324 


45675 


141 


< 0.1 


< 0.1 


1.4 


1.58518830 


28719 


57 


< 0.1 


< 0.1 


1.5 


1.67192622 


18964 


25 


< 0.1 


< 0.1 



rithm of the correlation length and thus remains finite 
It is therefore reasonable to think of using the half- 
chain von Neumann entropy, to detect a phase transition. 
It turns out that when running the iTEBD algorithm, S 
converges faster to a steady value than the mean value of 
the order parameter and thus provides a faster detection 
of the position of the critical point (varying the mag- 
netic field, A, and scanning for the peak of S). Yet, the 
von Neumann entropy converges more slowly than the 
energy, see Fig. [T31 Around the critical point, the spec- 
trum of the Hamiltonian is filled with a lot of low-energy 
excited levels which energy is very close to that of the 
ground state. An arbitrary superposition of such excited 
states will have energy close to that of the ground state, 
but can in principle exhibit very different entanglement 
properties. We believe that this is why it takes much 
longer to get a reliable estimate of the entropy. One has 
to make the energy converge close enough to that of the 
ground state so that the entropy of the obtained state 
also faithfully reflects that of the ground state. 



Energy^ | o Entropy | 



TABLE IV: Errors in the energy in relation to the exact value 
for X = 2,4,8 and 16. These errors are greater around the 
critical point (which for x = 2 is close to A = 1.1). These 
values were obtained with e = 0.1, showing a clear dominance 
of truncation error over the errors introduced by finite step 
time evolution. All entries of this table should be multiplied 
by 10-«. 



200 300 
imaginary Time 



A 


Exact 


X = 2 


X = 4 


X = 8 


X = 16 


0.5 


421292 


-3 


-3 


< 0.1 


< 0.1 


0.6 


914778 


-36 


-36 


< 0.1 


< 0.1 


0.7 


1869961 


-389 


-389 


< 0.1 


< 0.1 


0.8 


3804448 


-4255 


-4255 


-1 


< 0.1 


0.9 


8484551 


-66920 


-66920 


-12 


-2 


1.1 


47444179 


-437545 


-437545 


-6473 


-3 


1.2 


36551466 


-74387 


-74387 


-270 


5 


1.3 


30064632 


-20221 


-20221 


-23 


5 


1.4 


25539496 


-6976 


-6976 


< 0.1 


5 


1.5 


22144107 


-2797 


-2797 


4 


5 



TABLE V: Errors in the entropy for different values non- 
critical A and of x- All values have been multiplied by 10*. 
Note the increasing accuracy as a function of x- 

Let us now clarify why we were interested in reaching 
full convergence for the half-chain entropy. A common 
method to locate a phase transition is to analyze the vari- 
ation of an order parameter. On another hand, we know 
that the half-chain entropy of a critical system diverges 
while, off and close to criticality, it scales as the loga- 



FIG. 13: Convergence of the energy and entropy, at the effec- 
tive critical point, during the imaginary time evolution, with 
X = 8, A = 1.006 and e = 10^^. The full convergence of the 
energy (eight decimals) took ~ 10^ steps while ^ 6 10^ steps 
where necessary to make the entropy converge. 



APPENDIX B: METASTABILITIES 

An important issue when running the iTEBD algorithm 
is to be sure that one is not driven to a local minimum. 
Here we point out the existence of some meta stabilities 
in the simulation with respect to the choice of the initial 
state (an effect which is also present in standard DMRG 
simulations). In all our calculations, we have used an ini- 
tial state which matrices Ta and F b (see [l^l for details) 
are of the following form : random entries in the 2x2 
left upper corner, and all other entries set to zero. How- 
ever, when performing a simulation for the Ising chain 
for some value A of the transverse field, one could use, 
as initial state, the result of a simulation performed at 
some close value of A. Although this procedure can sub- 
stantially decrease the time necessary to make the energy 
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and the half-chain entropy converge, it can also lead to 
misleading results regarding the position of A* , taken as 
the point where the order parameter vanishes, as can be 
seen on Fig. 1141 Simulations which start from a previous 
minimization run of a larger A do produce unphysical 
results. Thus, all simulations must start from random 



I 0.3 

I 

0.2 



o Random 
□ Previous 
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FIG. 14: Magnetization as a function of the transverse mag- 
netic field using different initial states for x = 4 ''■nd e = 0.1. 
In one case (open circles) we use a random 2x2 matrix as 
an initial state. In the second case (open squares) the initial 
state for A = Ao is the final state obtained for A = Ao -1-0.001. 
The two methods do not give similar results for the position 
of the critical point. 



APPENDIX C: BOOSTED ITEBD 

The performance of the iTEBD algorithm depends on 
the initial conditions and the gap above the ground state. 
The results of our study suggest that using finite x one 
is perturbing the system in a way similar to have an ef- 
fective gapped Hamiltonian. However, if the gap is small 
the convergence of the algorithm can be very slow. To 
see this, we can consider as initial state a state with non 
zero projection on the ground state 



(CI) 



with |a| < 1. It is easy to see that, if the Hamiltonian 
has a gap A, the Euclidean evolution of an initial state 
with non zero projection on the ground state will lead to: 



IV') = exp(-i?T)|V) = aexp(-i?or)|V'o) + v/(l-"')lV'l) 

(C2) 



with ji/;^) — exp{—Ht)\ip±). From 
(VAI^IVA) A + i?o, 



we see the long time limit of the above expression, differs 
from the ground state (as already pointed out in ref. [H) 
by terms of the order: 



KV-olV')! 



1 



2a2 



■ exp(-2TA). 



(C4) 



Now if we approach the critical point of a phase transi- 
tion we know that the correlation length scales with the 
critical index v of the corresponding universality class 



(C5) 



(C3) 



where t denotes, again, the distance from the critical 
point. Assuming that A ~ 1/^, we see that, even in the 
case of a good guess of the initial state (that is, in the case 
^ 1), the convergence of the algorithm slows down in 
the critical region. In order to partially cure this slowing 
down we can perform a linear extrapolation of the results 
obtained after a small interval of Euclidean time dr and 
get a new estimate for the MPS. Given a generic element 
of the MPS matrix at Euclidean time r, A{si){T), and 
the same element at time t -|- dr, A{si)(T + dr), we con- 
struct a new MPS which matrix elements A{si){T) are 
the extrapolation at time T of the straight line passing 
from the two points at r and t + dr. Before promoting 
the guessed MPS to a new initial condition, we should 
check that the state it describes has a lower energy than 
the MPS obtained a,t t + dr before the extrapolation. 
A lower energy indeed means a greater overlap with the 
ground state. In this case the extrapolation is successful 
and we promote the guess to an initial condition of the 
new evolution. In case the energy of the new extrapo- 
lated state is greater than the energy of the state before 
the extrapolation, we just neglect it and keep the state 
we had before the extrapolation. The new Euclidean evo- 
lution is also of length dr and is followed by the attempt 
of a new extrapolation. We iterate the procedure till 
we reach convergence. We call this technique boosted 
iTEBD. 

In order for the extrapolation to work, we have to tune 
finely its two parameters: the waiting time dr and the 
amount of time we extrapolate, T. Once these two pa- 
rameters are fixed, we are able to accelerate the conver- 
gence by a factor greater than 10. A typical case is shown 
in Fig. [15] where we show the convergence of the simula- 
tions of a X = 32 MPS at A 1. Without the boosted 
iTEBD algorithm, with e = 0.01, after the number of Trot- 
ter steps considered, the system had still not converged. 
Increasing e = 0.1 translates in a coarser precision but 
with a convergence time about ten times shorter. A fur- 
ther improvement in convergence is obtained by keeping 
the same e — 0.01, and hence the same precision, but 
boosting the evolution with the extrapolation technique 
described. We see that the gain in convergence time is 
bigger by a factor of ten. As we can see, there is a point 
where the extrapolation fails and the normal evolution 
is continued. We can also check that before the first 
extrapolation, the boosted evolution coincides with the 
unboosted evolution with the same e. 
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FIG. 15: Boosted itedb algorithm. We plot the half-chain 
entropy as a function of the Trotter steps for a x = 32 MPS 
at A = 1. This is taken as a typical case from a large number 
of examples with different x and magnetic fields that present 
a similar behavior. We compare the results obtained with 
£ = 0.01 with both boosted and standard ITEBD and the re- 
sults obtained with e = 0.1 with standard iTEBD. As we can 
see, the unboosted case with e = 0.01 is far from having con- 
verged in the number of Trotter steps considered. On the 
other hand, the boosted system with e — 0.01 converges (up 
to 8 decimals) in a smaller amount of time steps than the un- 
boosted algorithm with an e ten times bigger. Indeed, the lat- 
ter simulation has still not converged in the window of Trotter 
steps shown. The discrepancy in the asymptotic values is due 
to e corrections described in the previous appendices. From 
this plot, we can safely deduce that the effect of the boost is 
to reduce the convergence time by a factor greater than 10 in 
the case we have analyzed. 



APPENDIX D: COMPARISON WITH DMRG 

In order to ensure that the effects we observe are not 
artifacts of the algorithm used, we reproduced some of 



them with a different algorithm. We have chosen to use 
the open source code for DMRG written by the Pisa 
group (39| . This program performs an infinite DMRG 
update of the system by growing it till it reaches a cho- 
sen chain length. At this stage, it performs several finite 
size sweeps through the chain (at least three in our case) 
in order to compute the reduced density matrix of all pos- 
sible chain bipartitions and improve the infinite results 




Energy Entropy 
DMRG N=16284 -1.27321717 0.68557374 
iTEBD -1.27323939 0.68065196 

iTEBD - DMRG(N=16384) -0.00002222 -0.00492178 



TABLE VI: Comparison of DMRG energy and entropy results 
with the infinite size result produced by iTEBD with e = 
10~* and where both methods used x = 16. 



We checked the stability of the presented results on 
the variation of the number of finite size sweeps. In this 
way, we are sure that the results have converged. We 
have checked that for a fixed number of level (to in the 
language of DMRG that corresponds to x in this pa- 
per), increasing the chain length makes the results con- 
verge to those obtained with the algorithm we have used 
in the paper. It was interesting to see that the DMRG 
convergence is however quite slow as compared with the 
boosted iTEBD. 
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